##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(ggrepel)
library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)
library(ggsci)
library(ggrastr)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    gene <- "MUC6"

    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic"
    input_file <- paste0("~/20220915_gastric_multiple/dna_combinePublic/finalPlot/MUC6_BMP6_CFTR/Diff/DiffGene.Pit.tsv")
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6/")

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

out_path <- opt$out_path
gene <- opt$gene
input_file <- opt$input_file

dir.create(out_path , recursive = T)

##########################################################################################

data <- data.frame(fread(input_file))
data2 <- data
colnames(data)[3] <- "-log10_P_Value"

##########################################################################################

#设置阈值
q_t <- 0.05
foldchange_t <- 1
logFC_cutoff <- log2(foldchange_t)
log10_P_Value_cutoff <- -log10(q_t)

## 展示想展示的基因
data$gene <- ifelse( data$gene %in% c("GKN1" , "GKN2" , "DSC2" , "DSG2") , data$gene , "" )

col <- c('#376B6D','#9F353A' , '#F7C242','#BDC0BA' )
names(col) <- c('DOWN', 'UP', 'TARGET', "GREY" )

alpha <- 1

options(ggrepel.max.overlaps = 200)
plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
        geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["UP"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                aes(size = abs(log_FC)),col = col["DOWN"],alpha = alpha)+
        geom_point(data = subset(data,gene!=""),
                aes(size = abs(log_FC)),col = col["TARGET"],alpha = alpha)+
        theme_bw()+
        labs(x='log fold change\n(Scissor+ cells/All other cells)',y='-log10(adjusted p-value)')+
        geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
        geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
        geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff), 
            fontface="bold", nudge_x = 0.65, nudge_y = 1, 
            aes(label = gene),size = 5.5,col = 'black' , face = 'bold') +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='none',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 10,color="black",face='bold'),
                axis.title.x = element_text(size = 14,color="black",face='bold'),
                axis.title.y = element_text(size = 14,color="black",face='bold'),
                strip.text.x = element_text(size = 12,color="black",face='bold'),
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5)
                )

image_name <- paste0( out_path , "/DiffGene.Pit.pdf" )     
ggsave( image_name , plot1 , width = 5.5 , height = 6 )

##########################################################################################
## 通路富集

computGO <- function(gene_symbol = gene_symbol , type = type){
        gene_df <- bitr(gene_symbol,fromType = 'SYMBOL',toType = 'ENTREZID', 
                    OrgDb = 'org.Hs.eg.db')
        ego <- data.frame(enrichGO(gene_df$ENTREZID,OrgDb = org.Hs.eg.db,ont = 'ALL',readable =T))
        ego$type <- type
        return(ego)
}

result_diff <- data2
result_diff$CLUSTER <- ifelse(result_diff$log_FC > log2(foldchange_t) & result_diff$X.log10_P_Value > log10_P_Value_cutoff , paste0( "MUC6_Mut_high_expression"),
ifelse(result_diff$log_FC < -log2(foldchange_t) & result_diff$X.log10_P_Value > log10_P_Value_cutoff , paste0("MUC6_Mut_low_expression"),"no_sig"))
result_diff <- subset(result_diff,CLUSTER %in% c(paste0( "MUC6_Mut_high_expression"),paste0( "MUC6_Mut_low_expression")))

enrich_df <- c()
for( clust in unique(result_diff$CLUSTER) ){
        print(clust)
        gene_symbol <- unique(subset( result_diff , CLUSTER == clust )$gene)
        tmp_go <- computGO(gene_symbol = gene_symbol , type = clust)
        tmp_go$class1_vs_class2 <- "MUC6_Mut"
        enrich_df <- rbind( enrich_df , tmp_go )
}

write.table(enrich_df, paste0(out_path , "/enrich.go_filt.txt"), sep = '\t', row.names = FALSE, quote = FALSE)
enrich_df <- subset( enrich_df , 
        !Description %in% c(
                "cardiac conduction" , 
                "bundle of His cell to Purkinje myocyte communication" , 
                "regulation of ventricular cardiac muscle cell action potentia") 
        )

##########################################################################################
## 每种类型的通路取前5最显著的
result <- c()
top_n <- 5
for( class_type in unique(enrich_df$class1_vs_class2) ){
    for(t_type in unique(enrich_df$type)){
      for( t_path in unique(enrich_df$ONTOLOGY)){
        tmp <- subset( enrich_df , type == t_type & ONTOLOGY == t_path & class1_vs_class2 == class_type )
        tmp <- tmp[order(tmp$qvalue , decreasing=F),][1:top_n,]
        result <- rbind(result , tmp)
      }
    }
}


## 画图
plotBar <- function(result = result , typeN = typeN){
        result_tmp <- subset( result , type == typeN )
        result_tmp$log10Pvalue <- -log10(result_tmp$pvalue)

        p <- ggplot(result_tmp) +
                labs(title = typeN) +
                facet_grid(.~ONTOLOGY , scales = "free_x") +
                aes(x = Description, y = log10Pvalue, fill = ONTOLOGY) +
                geom_bar(stat = "identity",colour="black") +
                scale_fill_npg()+
                ylab("-log10Pvalue") +
                theme(
                axis.title=element_text(size=15,face="plain",color="black"),
                axis.text = element_text(size=12,face="plain",color="black"),
                axis.text.x = element_text(angle = 90,colour = "black",hjust=1,vjust=0.6),
                axis.title.x = element_blank(),
                legend.postion = "none" ,
                legend.background = element_blank(),
                panel.background = element_rect(fill = "transparent",colour = "black"),
                plot.background = element_blank()
        )

        return(p)
}

for( class_type in unique(enrich_df$class1_vs_class2) ){
        result_tmp <- subset( result , class1_vs_class2 == class_type )
        typeN <- unique(result_tmp$type)[1]
        p1 <- plotBar(result = result_tmp , typeN = typeN)
        typeN <- unique(result_tmp$type)[2]
        p2 <- plotBar(result = result_tmp , typeN = typeN)

        p <- p1 + p2
        out_file <- paste0( out_path , "/diffexp." , class_type , ".gobp." , foldchange_t , ".pdf" )
        ggsave(filename = out_file, plot = p, width = 15, height = 10)
}